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In this paper we study the use of cross-correlations between multiple gravitational wave (GW) 
data streams for detecting long-lived periodic signals. Cross-correlation searches between data from 
multiple detectors have traditionally been used to search for stochastic GW signals, but recently 
they have also been used in directed searches for periodic GWs. Here we further adapt the cross- 
correlation statistic for periodic GW searches by taking into account both the non-stationarity and 
the long term-phase coherence of the signal. We study the statistical properties and sensitivity of this 
search, its relation to existing periodic wave searches, and describe the precise way in which the cross- 
correlation statistic interpolates between semi-coherent and fully-coherent methods. Depending on 
the maximum duration over we wish to preserve phase coherence, the cross-correlation statistic can 
be tuned to go from a standard cross-correlation statistic using data from distinct detectors, to the 
semi-coherent time-frequency methods with increasing coherent time baselines, and all the way to a 
full coherent search. This leads to a unified framework for studying periodic wave searches and can 
be used to make informed trade-offs between computational cost, sensitivity, and robustness against 
signal uncertainties. 



I. INTRODUCTION 

Long lived quasi-periodic gravitational waves (GWs) 
from rapidly rotating non-axisymmetric neutron stars 
are among the promising sources of detectable GWs for 
ground based detectors such as LIGO, Virgo, GEO600 
etc. A number of searches for long-lived periodic GWs 
have been carried out using data from ground based GW 
detectors. These include searches using data from the 
interferometric and bar detectors. These searches are of 
two kinds depending on the size of the parameter space 
that is searched: 

i. Targeted searches for sources whose parameters are 

well known from other astrophysical observations 
[H m [3]. Such searches are not computationally 
intensive, and use statistically optimal matched fil- 
tering techniques. 

ii. Wide parameter space searches either for neutron 

stars in binary systems whose parameters are 
poorly constrained from prior observations [4l, or 
blind searches for as yet unknown neutron stars 

While none of the above searches have yet resulted in a 
detection, there have been some notable successes. For 
the searches targeting known pulsars, the limits on the 
gravitational wave emission and the corresponding limits 
on the deformation are starting to become astrophysi- 
cally interesting. 

Similarly, a lot of the groundwork has been laid for 
meeting the computational challenges for the wide pa- 
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rameter space searches. Computationally efficient meth- 
ods and hierarchical data analysis pipelines have been 
developed which allow us to vastly improve the ratio of 
sensitivity to computational cost. Most of these are semi- 
coherent methods, i.e. combinations of coherent analy- 
ses combined together by excess power techniques, and 
they come in two main flavors. The first combines short 
segments of simple Fourier transformed data. The base- 
line of the short Fourier transforms is chosen such that 
the signal manifests itself as excess power in a single fre- 
quency bin, and the excess power is combined by var- 
ious methods. The simplest is the StackSlide method 
[S] which adds the normalized excess power from the 
short segments, taking care to "slide" the frequency bins 
to account for the Doppler shift and intrinsic spindown. 
The PowerFlux method |T(J) is very similar; it performs 
a weighted sum of the normalized power using weights 
which take the sky-position and polarization dependent 
sensitivity of the detector into account; the weights serve 
to improve the sensitivity. Finally, there is the Hough 
transform method which performs a weighted sum of 
binary-number counts calculated by setting a threshold 
on the normalized excess power. This is more robust 
and computationally efficient, though at the cost of be- 
ing somewhat less sensitive. All three methods have been 
used to analyze LIGO data in all-sky wide frequency 
band searches for GWs from isolated neutron stars ^ ^, 
and these are so far the most sensitive wide parameter 
space GW searches of their kind published so far; we shall 
refer to them as the "standard" semi-coherent searches 
in the rest of this paper. 

A variant of these standard semi-coherent techniques 
are the so-called hierarchical searches which aim to 
search deeper by increasing the coherent time baseline 
[HI [HI [12] • This requires a sky-position (and spindown) 
dependent demodulation to be performed before calcu- 
lating the excess power statistic. The extra demodula- 
tion step significantly increases the computational cost 
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and such a search pipeline is currently being employed on 
larger computational platforms such as EinsteinOHome^. 

In addition to the above surveys for isolated neutron 
stars, searches have also been carried out for gravita- 
tional waves from neutron stars in binary systems. A 
plausible argument for why some neutron stars may be 
emitting detectable GWs applies to neutron stars in bi- 
nary systems, and in particular, to the Low Mass X- 
ray Binaries (LMXBs) which consist of a neutron star 
and a low mass main-sequence star. The observed X-ray 
flux from these systems is due to the high rates of accre- 
tion of matter onto the neutron star. It is observed that 
the rotation rates of neutron stars in LMXBs is signifi- 
cantly lower than that might be expected on theoretical 
grounds; the highest theoretically possible rotation rate 
is significantly larger than that of a IkHz, while the cur- 
rent observed record is ^ 620 Hz. It was suggested (first 
by Bildsten [IT) that this apparent upper bound on the 
rotation rate might be due to a balance between the spin- 
up due to accretion and the spindown due to the emis- 
sion of gravitational radiation - there is virtually a "wall" 
created by the flux of GW radiated, which increases as 
f2^, where fl is the angular rotational frequency of the 
spinning neutron star and this limits its spin-up. There 
are a number of other suggested explanations which do 
not involve gravitational radiation, but accreting neutron 
stars are clearly promising sources of detectable gravita- 
tional radiation. So far two searches have targeted Sco 
X-1, the brightest LMXB. These have used very different 
techniques; [4j used a coherent integration on 6 hours of 
data from the second science run of the LIGO detectors, 
while [7] uses a cross-correlation statistic on data from 
the more recent fourth science run. The elucidation and 
generalization of this cross-correlation technique tailored 
to periodic GW searches, and its relation with the other 
searches discussed above will occupy us for the rest of 
this paper. 

The results from these searches are starting to become 
astrophysically interesting. For example, using data from 
the latest science runs of the LIGO detectors, it is ex- 
pected that the indirect spindown limit on the ampli- 
tude of gravitational waves from the Crab pulsar will be 
beaten by about a factor of 3. The resulting limits on the 
ellipticity of the known pulsars are also starting to place 
constraints on the equations of state of nuclear matter in 
neutron stars (see e.g. [T31[TS]). A detection would lead 
to new insights about neutron star physics not obtain- 
able by other means. Searches using large amounts of 
data from the LIGO detectors operating at design sensi- 
tivity are well underway, and the results are expected to 
become yet more astrophysically interesting in the near 
future. 

Almost all of these searches mentioned above have been 
based on techniques which look for signals of a given form 
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in a single data stream, i.e. either matched filtering tech- 
niques or semi-coherent power summing methods. While 
both matched filtering and semi-coherent techniques have 
been generalized and used to analyze data from multiple 
interferometers [HI \W\ , the starting point for these meth- 
ods is always the analysis of a single data stream. There 
is however one exception, which is the method used in 
[3 I17| and is inherently based on looking at multiple 
data streams. Let us consider two data segments 

xi{t) for te[Ti- AT/2, Ti + AT/2] , (1.1a) 
X2{t) for t e[T2- AT/2, T2 + AT/2] . (1.1b) 

If a signal resulting from the same gravitational wave is 
present in both streams, it should be possible to cross- 
correlate the output of two detectors to extract the signal. 
The basic cross-correlation statistic is 

<.Ti+AT/2 1.T2+AT/2 

/ dti / dt2Xi{ti)x2{t2)Q{tl,t2) , 

JTi-AT/2 ■/T2-AT/2 

(1-2) 

where Q(ti,t2) is an appropriately chosen filter function. 
This technique was originally developed for the stochastic 
background searches where the cross-correlation is abso- 
lutely essential and is based on the fact that multiple 
detectors will see the same GW signal [HI [19], and it 
has been used extensively to search for a stochastic GW 
background using LIGO data [SnUUKS^. The function 
Q{ti,t2) can be tuned to search for GWs coming from 
a particular sky position and also polarization [19] and 
this method has been used to search for periodic waves 
from the neutron star in Sco X-1. All previous discus- 
sions of this method have however been in the context 
of stochastic searches. In this paper, we investigate in 
detail its applications for periodic wave searches. 

The optimal form of the function Q(ti, t2) depends on 
the kinds of sources that we are looking for. Thus for a 
stochastic background we use the facts that the statistical 
properties of the signal are time independent and that 
the two polarizations are statistically independent. In 
particular, the optimal Q is time invariant, i.e. a function 
of only the difference ti — ^2- Furthermore, Q turns out 
to depend on the expected spectrum of the stochastic 
background. 

For periodic GWs from neutron stars, many of these 
assumptions do not hold. The signal is deterministic and 
non-stationary (because of the Doppler shift), and the 
two polarizations are not independent. There is yet an- 
other ingredient present for periodic signals that is not 
present for stochastic sources. In principle, since the sig- 
nals we are looking for have long term phase coherence, it 
should be possible to cross-correlate any pair of data seg- 
ments to extract the signal, regardless of how far apart 
the segments are in time and regardless of whether they 
are from the same interferometer or not. It will turn 
out that the sky-resolution is much coarser than for the 
standard periodic searches; the appropriate baseline is 
not the Earth-Sun distance but rather the distance be- 
tween the two detectors. This leads to a much lighter 
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computational burden for a blind search. All of these 
issues will be discussed in detail in the rest of this paper. 

The paper is organized as follows. Sec. |TT]sets up no- 
tation and describes the waveforms that we are looking 
for; this includes both isolated neutron stars and neu- 
tron stars in binary systems. It also discusses the short 
segment Fourier transforms (SFTs) and the restrictions 
on their time baseline for the signal power to be concen- 
trated in a single SFT frequency bin. Sec. Ill motivates 
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and defines the basic cross-correlation statistic for a pair 
of short data segments; Sec. |IV| discusses the statistical 
properties and the sensitivity of the search; Sec. |V] elu- 
cidates the relation of the cross-correlation method with 
the !F statistic; Sec. |VI] provides estimates of the parame- 
ter estimation that can be achieved and Sec. |VlT] investi- 
gates the question of resolution of parameters such as sky 
position, spin-down etc. Sec. VIII concludes with a sum- 



mary of our results and suggestions for future work, and 
finally appendix discusses some technical and concep- 
tual issues which have been ignored in the earlier sections 
for simplicity. 



II. NOTATION AND USEFUL EQUATIONS 

A. The waveform 

The waveform we are looking for is a tensor metric 
perturbation 



h{t) = h+{t)"e+ + hy,{t)"ey 



(2.1) 



where {eU|A = +, x} is a transverse-traceless polariza- 
tion basis associated with the GW propagation direction 
and tailored to the polarization state of the waves so that 

h+{t) ^ A+ cos <^{t) , /ix(t) = Ax sin$(t) . (2.2) 

If L is the angle between the line of sight n to the star 
and its rotation axis, the amplitudes are 



^4 



1 + COS^ i 



hoAy 



(2.3a) 
(2.3b) 



In the neutron star rest frame with proper time t, the 
phase is 



<i>(t(T))=$0 + 2^<i/0T+^/lT2. 



(2.4) 



The reference time where all the spindown parameters 
are defined is taken to be r = 0, and <i>o is the phase at 
T = 0. 

A detector's scalar strain response is the contraction 
of the tensor metric perturbation with a response tensor^ 



h{t) = h{t) -.lit) = J2 FA{t)hA{t) 



A=+,x 



where 



FA{t)^"eA ■■ d{t) 



(2.5) 



(2.6) 



The polarization basis {j^a} is sometimes inconvenient, 
because its definition involves not only the direction to 
the source but also the source's polarization state (specif- 
ically the orientation of the neutron star's spin). For a 
given sky direction n, one can always construct a trans- 
verse, traceless polarization basis by starting e.g., 
with the vector transverse to ft and lying in the Earth's 
equatorial plane. The relationship between this reference 
basis and the preferred polarization basis of the source is 
described by the polarization angle ?/;: 



= cos 2'i/; + sin 2'!/' 
Vx = — e+ sin 2-0 + e'x cos 2ip 



That means, if we define 



a{t;n) = d{t) : £>(«) 
b{t; n) = d{t) : *ex (n) 



(2.7a) 
(2.7b) 



(2.8a) 
(2.8b) 



(which are time-dependent because of the rotation of the 
detector tensor d), we can decompose the beam pattern 
functions as 

F+{t;n,tp) ^ a{t;n) cos 2iIj + b{t;n) sm2tp , (2.9a) 
i^x {t; ft, Tp) = b{t; n) cos 2-0 - a{t; n) sin 20 . (2.9b) 

The polarization angle is a property of the source, but 
the functions a{t; n) and b{t; n) depend on both the sky 
position of the source and the detector in question. 



1. Isolated neutron stars 

The relation between the detector time t and the neu- 
tron star time r depends on whether the neutron star is 
isolated or in a binary. For an isolated neutron star, we 
assume^ that the star is at rest with respect to the SSB 
frame. Let f{t) be the position of the detector in the SSB 
frame and v{t) its velocity. The times of arrival of the 
wave at the detector and the SSB are 



t 



relativistic corrections 



(2.10) 



^ For an interferometer with arms along the unit vectors u and v, 
d = ^{u (S) u — V (S) v). 



As it turns out, so long as the neutron star is moving inertially, 
this assumption is not necessary; the frequencies involved are all 
simply offset by the constant Doppler shift between the neutron 
star rest frame and the SSB. 
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The relativistic corrections can be ignored for our pur- 
poses. The instantaneous frequency is then, to a very 
good approximation 



Roemer delay: 



i = T- 



-relativistic corrections . (2.14) 



fit) = m + fit)- 

f{t) = h + ht. 



(2.11) 
(2.12) 



The parameters of the signal from an isolated neutron 
star are thus the so-called amplitude parameters (or 
nuisance parameters) {/lo, cos i, -0, <i>o} and the Doppler 
parameters A = {n, /q, /i, . . .}. The Doppler param- 
eters determine the frequency evolution of the signal 
through ( |2.11[ ). The frequency and spindown ranges will 
canonically be taken to be 50 Hz < /o < 1000 Hz, and 
— 1 X 10~*Hz/s < /i < 0. These were the ranges used 
in [B]. The lowest frequency is determined by the per- 
formance of the detector, and it will be lower for the ad- 
vanced detectors. The upper end of the frequency range 
could conceivably be as high as 2000 Hz depending on the 
computational cost. 

Written in terms of the detector time and including 
first spindowns, the phase is: 



1 



r • n 



$(t)=$o + 2^ (^/oi+^/it'j +27r(/o + /ii)- ^ 

(2.13) 

We have ignored the \ fi{r- njcf' term. In fact, even the 
term fit{r-n/c) will be ignored in most of the calculations 
below.'' 

Let us quantify the restrictions on the parameter space 
due to these approximations adapting the "1/4-cycle 
criterion" used in |52j: any physical effect which con- 
tributes less than 1/4 of a cycle to the phase of the sig- 
nal over a given coherent observation time will be ig- 
nored. Since \r-n/c\ < lAU/c w 500s, we will have 
i|/i|(f • r?/c)2 < 1/4 if l/il < 2 X 10-6 Hz/s. This 
is much larger than any spindowns we can realistically 
consider. On the other hand, the fit(f- n/c) is, in gen- 
eral not negligible for realistic spindowns and observation 
times of months. However, we will break up our obser- 
vation time into shorter segments of duration much less 
than a day. Over say 1 hour, this term is ignorable if 
|/i| < 3 X 10~^Hz/s which is still a very large spindown. 



2. Neutron stars in binary systems 

To account for the motion of the neutron star in a 
binary orbit, we need to add the orbital time delays to 



(2.101. The most important contribution is again the 
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approximations, and nor do they ignore the relativistic Einstein 
and Shapiro corrections. 



Here Forb is the position vector of the neutron star in the 
binary system's center of mass frame. 

There are four relevant orbital parameters. The first 
is the orbital period Porb, and, if available, its derivative 
-Porb- We then need a reference time within the orbit for 
which we use Tasc , the time of crossing of the ascending 
node. The third parameter is the projected semi-major 
axis of the neutron star, Op = a^sini. The final param- 
eter is the orbital eccentricity e. In addition, there are 2 
parameters specifying the orientation of the orbital plane, 
i.e. the inclination angle i (not to be confused with the 
orientation of the neutron star axis l) and the argument 
of periapsis lu. Of these 6 parameters, only 5 are required 
to define the phase model because of the projection along 
the line of sight n; see p4] for further details. 

We therefore have a total of 5 parameters of the bi- 
nary which determine the frequency evolution of the sig- 
nal: Abin = (dx sin i, e, Porb, ^asc, w). In the case when 
the orbit is circular (e = 0), the argument of periapsis 
and the initial orbital phase combine additively into a 
single parameter so that we are left with only 3 search 
parameters: Abin = (op, Porb, Tasc)- We will not include 
higher derivatives of Porb- As an example, for Sco X- 
1 (the brightest LMXB), some of the orbital parameters 
are Porb « 6.8x lO'^s, and Up/c « 1.44 s, and e < 3x 10"^ 



Let i/orb be the velocity of the neutron star in the 
center-of-mass frame of the binary. The observed fre- 
quency is, to a very good approximation, given again by 
the non-relativistic expression. 



m = fit) + m 



(2.15) 



Since v^rh is usually much larger than the Earth's orbital 
velocity, Worb is the dominant contribution to the Doppler 
shift. 



B. Short-time Fourier transforms 

Given a time series detector output from a detector, 
it is convenient to break it up into short segments of 
length AT and to store the Short-time Fourier Trans- 
forms (SFTs). The value of AT is chosen such that the 



approximation (2.24 1 is valid and as we will see, this leads 



to different restrictions on AT for neutron stars which are 
isolated or in binary systems. Such SFT databases are 
commonly used in the LIGO, GEO and Virgo collabo- 
rations for periodic wave searches, and we will also base 
our data analysis strategies mostly on SFTs [27] . 

Let x{t) be a time series sampled discretely at intervals 
of St. Let us consider N samples Xj for j = . . . — 
1, and let AT = N5t. Our convention for the discrete 
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Fourier transform will be 

JV-l 

Xk = St Xje' 

3=0 



i2TTjk/N 



(2.16) 



where fc = 0, 1 . . . (TV - 1). For < fc < [N/2\, the 
frequency index k corresponds to a physical frequency 
fk = k/AT with [.J denoting the integer part of a given 
real number. The values lN/2\ < k < N — 1 correspond 
to negative frequencies given by fk — {k — N)/AT. Each 
SFT stores the real and imaginary values of Xk for a range 
of frequency bin indices k. The SFT will span the 
time interval [Tj - AT/2, Ti + AT/2]. When necessary, 
we will denote the data at the fc*^ frequency bin of the 
SFT by ikj- 

Eq.(2.16) is actually a simplification. In practice, to 



avoid spectral leakage, a taper wj is applied while taking 
the Fourier transform: 



3=0 



^-i2TTjk/N 



(2.17) 



See e.g. [2S] for details. We will mostly ignore window- 
related issues in this paper. 

The detector output x{t) is the sum of noise n{t) plus 
a possible gravitational wave signal: 



x{t) = n(<) + h{t) . 



(2.18) 



C. The short-duration Fourier transform of the 
signal 

We now calculate the Fourier transform of the signal 
over an observation duration [T — AT/2,T + AT/2] cen- 
tered at the time T. We assume AT is small enough so 
that {-FaI^ = +, x} can be treated as constants in this 
duration; this means AT <C 1 day. We assume that the 
observation duration is small enough so that the phase 
of the signal in this duration can be expanded in a power 
series at the mid-point T: 



$(t) = $(T) -I- 27r/(r)(t - T) . 



(2.22) 



The validity of this approximation sets the limits on how 
large AT can be. If /(<) is the time-derivative of the 
signal frequency at any given time t, the above approxi- 
mation is valid whenever effects of the frequency deriva- 
tive / can be ignored over the duration AT. Using the 
1/4-cycle criterion, this leads to / < AT~^. 

For isolated neutron stars, the time variation of f{t) 
is given by (2.111 and is due to two effects: the intrinsic 



spindown of the star, and the Dopplcr modulation due to 
the Earth's motion. Consider first the intrinsic spindown 
fi. Taking the largest spindown to be 10"*Hz/s, we get 
AT < lO^'s. For the Dopplcr shift, we can estimate / by 



keeping / fixed and differentiating vin (2.11 1. The result 



is worked out in [T^] and yields the following restriction 
on AT: 



We will assume the noise to be a real stochastic process of 
zero mean, stationary and Gaussian; in practice, we only 
need stationarity over a period AT, the time baseline 
of the SFTs. The properties of the noise are thus fully 
described by a single-sided power spectral density Sn{.f) 
which, in the continuous time case is defined as, 

/oo 
{n{t' + t)n{t'))e-'^^f*dt , (2.19) 
-oo 

where (•) denotes an average over an ensemble of noise 
realizations. Note that the average {n{t' + t)n{t')) is in- 
dependent of t' because of the assumption of stationarity. 
In practice, we are of course only given x{t) and not n{t) 
itself. So we must take care to ensure that the estima- 
tion of Sn{f) is not biased by the presence of a signal. 
Finally, the following expression for Sn is useful: 



(life 



^ |2\ 



AT 



Snifk) 



(2.20) 



This equation relates the variance of the (real and imagi- 
nary parts) of fik to the PSD, thus providing a more intu- 
itive understanding of the PSD. This is a special case of 
a more general expression which, in the continuous case, 
reads. 



1 



(2.21) 



AT < 4 X 10^ s X 



/ 500 Hz 
fo 



(2.23) 



In this paper, for isolated neutron stars, we will mostly 
use AT — 30 min as a canonical reference value. This 
is well within the above restrictions. The limits on AT 
are far more stringent for neutron stars in binary systems 
because of the higher Doppler shifts. The Sco X-1 search 
in m used AT = 60 s. 



With the approximation (2.221, in the time interval 



[T - AT/2, T + AT/2] we have 

h{t) =F+A+ cos($(T) + 27r/(T)(i - T)) 

+ Fx Ax sin($(T) + 2^/(T)(i - T)) . 

The Fourier transform of h{t) is easily seen to be, 

j.T+AT/2 

h{f) = / /l(t)e-^2./(t-T+AT/2)^^ 
JT-AT/2 



(2.24) 



e'-^(^) ^^+^+:^^-^-^ ^AT(/-/(r)) 



(2.25) 



where we have defined the finite time approximation 
'5at(/) := sin(7r/AT)/7r/ to the delta function 5{f). 
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This definition of the function 5at(/) leads to signifi- 
cant spectral leakage of the signal power into neighboring 
frequency bins. This can be improved by using suitable 
tapers as in (2.171. We assume that this has been done 



and we will henceforth assume that spectral leakage is 
neghgible. 



III. THE CROSS-CORRELATION STATISTIC 
FOR A PAIR OF SFTS 

Let us assume that we have two data streams cover- 
ing the time intervals X/ and Ij centered on the times 
T/ and Tj respectively; both intervals have the same du- 
ration AT. The data streams in the two intervals xi 
and xj could come from the same or different detectors, 
though of course if Tj = Tj then the detectors have to 
be different. The received signals in the two intervals are 
denoted by hi{t) {t e X/) and hj{t) {t € Xj) respectively. 
As before, we assume that the duration AT of the time 
intervals is such that the beam pattern functions are ap- 
proximately constant. We denote the PSDs of the noise 
in the two intervals by Sn\f) and Sn'^\f) respectively. 

The basic cross-correlation statistic corresponding to a 
filter function Q is, 

l-Ti+AT/2 i-Tj+AT/2 

Sij= dt dt' xi{t)xj{t')Qij{t,t') . 

JT,-AT/2 JT/-AT/2 

(3.1) 

We would like to understand how the optimal Qu can be 
chosen. The optimal choice depends in fact on the kind 
of signals we are looking for. The analysis presented in 
T8] describes the optimal choice of Q for stochastic sig- 
nals, and here we will tailor our discussion to the periodic 
signals described earlier. 

To get some intuition on the nature of Sij , let us eval- 
uate Sij in the frequency domain assuming that Qu is 
time invariant: Q{t,t') — Q(t — t'). Keep in mind how- 
ever that this will not be the optimal solution, and a 
more detailed analysis will be presented later. 



It is easy to evaluate (3.1) by writing xi{t) in terms 



of its Fourier transform. Along the way we approximate 
<5at by the delta function, but we however should not 
take Qij{t) to be a rapidly decreasing function of r as in 
[18] . Since our signals have long term phase coherence, 
Qij{t) will also turn out to be periodic. In any case, we 
still end up with the simple expression, 



Sij = 



dfi*i{f)i.Af)QijU) 



(3.2) 



The mean value of Sjj over an ensemble of noise realiza- 
tions is, 

/oo 
dfh}{f)hj{f)Qjj{f) . (3.3) 

Here we have assumed that the noise has zero mean, and 
that ni and nj are uncorrelated. If we assume further 



that hi <^ ni then the standard deviation is approxi- 
mately: 



AT 



df S^!\f)Sli^\f)\QjjU)\' . (3.4) 



Furthermore, it can also be shown under the same as- 
sumptions, that Sij and Sjk are uncorrelated for K ^ I: 



(SijSjk) — SikCTij 



(3.5) 



Thus, the correlation pairs formed from all pairs of dis- 
tinct SFTs are statistically independent. Note however 
that the same is not true for the third order moments; 
for example {SijSjkSki) 7^ even when the small signal 
approximation is valid. This is however not a problem 
for us because we will never need to calculate the third 
and higher order correlations between the 

Eq.(3.3l clearly demonstrates that taking Qij{t,t') 
to be time-invariant is, in general, suboptimal for the 
data analysis problem at hand. The signal frequencies 
// = f{Ti) and fj — f{Tj) at the midpoints of the two 
intervals are given by (2.11[) (for an isolated system) or 



(2.151 (for a binary systems) . In general, // and /,/ may 



be quite distinct from each other, especially if the in- 
tervals are far apart in time. Our assumptions on AT 
ensure that the signal power to be concentrated mostly 
in a single SFT frequency bin. Thus, no matter what we 
choose for Qij{f), the overlap between hj and hj might 
be quite small. This will lead to a small fijj and thus 
a small signal-to-noise ratio ^iij /crij. The fix is obvious: 
we need to shift the frequencies while constructing the 
cross-correlation statistic. So, if we define 6Jij = fj — fi 
then. 



S 



IJ 



dfS:*i{f)ij{f + Sfij)Qij{f + dfij/2) . (3.6) 



In the time domain, this corresponds to the non-time 
invariant filter: 



Qij{t,t') = e-^<'f"^^'+''^Qij{t-t') 
The mean fijj becomes, 



(3.7) 



^^1J■■={SIJ)^ / dfh*jif)hjif+dfij)Qij{f+Sfij/2), 

J —00 

(3.8) 

and the variance ajj is unchanged. 

An important quantity for us is the signal cross- 
correlation h*j{f)hj{f + 6fjj). We extract the amplitude 
term /iq and the delta-functions to define for / > 0, 



h*i{f)hj{f + Sfij) = hlGijSlAf - fi) ■ 



(3.9) 



The signal cross-correlation function Gu is an impor- 
tant quantity, much like the overlap-reduction function 
for stochastic searches defined in [T8j (though Qu is not 
exactly analogous to the overlap reduction function). 
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Apart from the frequency / and Tj,Tj, Qjj is a func- 
tion of the signal parameters, i.e. the ampHtude parame- 
ters {/iQ, i, V'j *I*o}j the Doppler parameters A, and possi- 
bly the binary parameters Abin- To avoid clutter, we will 



often drop the dependence of Qu on the signal parame- 
ters and T/, Tj, and just write Gij. 

Using (2.25) it is easy to calculate Gij- For / > 0, the 



dominant contribution is: 



Gu 



1 



{{Fi+Fj+Al + Fi^Fj^Al) - i{Fi+Fj^ - Fi^Fj+)A+A,, } , (3.10) 

A^ij = ^j{Ti)~^j{Tj) . (3.11) 

Here we have added the subscript / and J to the phase $ to emphasize that <i> is detector dependent. For isolated 
neutron stars, with the approximations explained in Sec. II A 1 this leads to 

Afjj ■ n 



k=0 



(3.12) 



We have used (2.13), ignored the fit{f ■ n/c) term, and 
defined Af/j i^T^T/) - r{Tj). Recall from that 
X are the same as x but without the factor of /iQ. 
We can now also average over cos i using the following 
relations: 



(^+^x)cos,. = . 



(3.13a) 
(3.13b) 



The average of Gij over cos l is thus, 
1 



/J / cos I 



60 



"{lFl+Fj++bFly,Fjy,) . (3.14) 



We can easily perform another average over the polariza- 



tion angle tp using (|2.9a 



{FixF.Jx)ip 

diabP 



-{aiaj + hjhj) 



(3.15) 



cd d'j , 



where, 



P™fd-l E ^A'^'^^A-l E ' (3.16) 



A=+, 



A=+, 



is a projection onto symmetric traceless tensors trans- 
verse to ft. This leads to: 



{Gu) 



1 



cos L^xjj 



{aiaj + bibj) 



— — e 
10 

— p uiab 0-j r- e 



(3.17) 



In the case of time-coincident SFTs, since A$/j reduces 
^niiiL^ this is just a normalization factor times the over- 
lap reduction function which would be used for a search 
for a stochastic background coming from a single point 
on the sky. [3 nil 121 130] 



IV. STATISTICS AND SENSITIVITY 

For each SFT pair (labeled by an index pair /J, we 
define the raw cross-correlation as the complex random 
variable: 



yk. 



ij 



Ar2 



(4.1) 



The frequency bin k' is shifted from k by an amount cor- 
responding to Sfjj: k' — k + [ATSful ■ Note that ykjj 
is computed using only data from single frequency bins 
in the two SFTs; this works under the assumption that 
the signal power is mostly concentrated in a single fre- 
quency bin. We emphasize that, this is not a fundamental 
limitation because we could, if we wished, consider the 
(optimally weighted) power from the neighboring bins 
as well if necessary. In the rest of this paper, we shall 
consider AT sufficiently small so that this assumption is 
valid. See Sec.|TT]for quantitative estimates on AT. 

In this section we initially make two additional sim- 
plifying assumptions. First we take the signal to be 
much smaller than the noise, i.e. h <^ n, and second 
we only consider yk.ij for / ^ J.^ The results ob- 
tained using these assumptions are probably the most 
relevant for practical applications. Firstly, for the ground 
based detectors the signal is indeed expected to be much 
smaller than the noise. Secondly, the number of pairs 
of distinct SFTs is much more than the number of self 
pairs; there is thus no significant loss in sensitivity if the 
self-correlations are not considered in the final detection 
statistic. 

The {yk.ij} are random variables with mean and vari- 



Both of these assumptions will be relaxed in Appendix [a| 
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ance given by, 



^kJJ 



'k,I.J 



1 

' 4Ar2 



= hlGij , 



(4.2) 
(4.3) 



To derive the expression for the mean, we have replaced 
^Arif — fi) by i5at(0) — AT, and for the variance we 
have assumed that the real and imaginary parts of Xk are 
uncorrelated and have the same variance. The {3^fc,7j} 
are not Gaussian variables, but we will only need their 
mean and standard deviation. 

Where convenient, we will replace the pair IJ with a 

single lowercase Greek index a,/3 Thus, yu.ij will 

often be denoted yk,a- To avoid unnecessary clutter, 
we also avoid putting the frequency index k explicitly in 
yk.a- In any case, one expects the signal contribution to 
be limited essentially to a single frequency bin k. Our 
task is now to combine the in a statistically optimal 
way to extract the signal amplitude /iq. The following 
analysis is very similar to what is used in [B] (see also 

mm)- 

We consider detection statistics which are weighted 
sums of the J^q : 



(4.4) 



We are interested in the probability distribution of the 
random variable p because this is required for comput- 
ing the sensitivity at given false alarm and false dismissal 
rates. It is simply obtained by examining the behavior 
of the noise in (of which p is made up of) which is 
derived from (4.1 1 by replacing the data x by the noise 



n in each data segment /, J. If we assume that the noise 
in each detector is Gaussian with mean zero, the noise 
in /5 is a sum of products of real independent Gaussian 
variables each having mean zero. Although J^q, is com- 
plex, the statistic p is real. The product of two inde- 
pendent Gaussian variables whose mean is zero, is a ran- 
dom variable whose probability density function (PDF) 
is essentially Kq{x), where K^^x) is the modified Bessel 
function of the second kind of order zero - more specifi- 
cally, if X - 7V(0, ax) and Y - iV(0, cry), then the PDF 
oi Z = XY is Kq{\z\/ (Jx<yY) /T^f^xCTY ■ This distribution 
has zero mean and a finite variance, namely, a^Oy- Then 
a generalization of the central limit theorem states that 
the sum of a large number of such zero mean variables 
tends to a Gaussian random variable [35]. Thus p is a 
Gaussian random variable whose mean p, and variance 
are given by: 

P = ^{UaPa + = ^1 ^{UaGa + <^a) , (4.5) 

a a 

a^ = 2^K|V^ (4.6) 

a 

Let us set a threshold pth on p to select detection candi- 
dates based on a false alarm rate a. It is easy to show 



that for Gaussian noise the threshold must be: 

pth = \/2aerfc-i(2a) , (4.7) 

where erfc is the complementary error function. The de- 
tection rate in the presence of a signal is. 



7 = - erfc I ^ 
2 V V2ct 



(4.8) 



Since /i oc /iq, this can be inverted to give the smallest 
value of h[j that will cross the threshold at given false 
alarm and detection rates. 



hf, = 25 



(4.9) 



where S — erfc ^(2a) — erfc ^(27). This can also be 
written in terms of the false dismissal rate /3 = 1 — 7 as 
S = erfc"^(2a) -herfc"^(2/3).'^ The solution for u„ which 
minimizes Hq can then be shown to be*". 



Ua CX 



(4.10) 



It is shown in appendix [A] that this solution also holds 
when we include the self-correlations (still assuming h <C 
n). 



Substituting from (4.101 into (4.4 1, the optimal detec- 
tion statistic is: 



p oc 



E 



y^Gi + y*o,ga 



(4.11) 



Substituting the expression for Ua from (4.101 back into 



(4.9), the optimal sensitivity is seen to be, 

\ 1/4 



ho 



(4.12) 



In the case when we are correlating data from two distinct 

interferometers with stationary noise floors Si^\f) and 
(2) 

if) J then a a is independent of a and is given by. 



1 



4AT2 



sl^Hf)si?Hf) 



(4.13) 



" This is proved by using the following property of the com- 
plementary error function; erfc(— a;) = 2 — erfc(2:). Setting 
X = — erfc~-'(27), we get 2 — 2/3 = 27 = erfc(— = 2 crfc(i'), 
which yields x = erfc" 1(2/3). 
^ This is perhaps easiest to see if we define a positive-definite inner- 
product over vectors x = {xa} as x ■ y := Re [x* yd cr^. In 

Hull 



terms of this inner product (4.9 I can be written as ho = -S 



where Ha = Ga/'^a- ^0 is then minimum when u is parallel to 
H. 
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We are using the superscripts in Sn'' and S^^'' to refer to 
the two detectors. Similarly, if we denote the average of 
over pairs of SFTs by then, 



\Qaf — ^pairs(|Sa|^)a 



(4.14) 



where iVpairs is the total number of SFT pairs. This leads 
to, 



ho = 



51/2 



V2{\g.?)'J' Kit, 



A 



0(1)0(2) 

I On On 

at" 



1/2 



(4.15) 



Similarly, if we relax the requirement that the pairs have 
to be from the distinct detectors, and instead assume 
that the noise floor in all SFTs is the same, Sn, then 



ho 



51/2 



1 



a/4 ^1/4 

/" pairs 



Sn 

AT 



(4.16) 



These are the equation we were after. They give us the 
sensitivity of the cross-correlation search as a function 
of the statistical false alarm and false dismissal rates, 
the SFT basehne AT, the noise floors of the SFTs, the 
number of SFT pairs iVpairs, and the geometrical factors 
contained in Qa. They tells us that the sensitivity grows 
coherently with AT and incoherently with -/Vpairs- Note 
however that we can correlate any SFT pair we like, so 
that iVpairs can be made much larger than the number of 
SFTs Nsft (even if we were to exclude self-correlations) . 
In fact, if we believe the signal to maintain phase co- 
herence over the entire observation time (which may be 
months or years), and if we can afford to do so compu- 
tationally, then TVpairs ~ N^ft so that ho cx {NshAT)''^/'^ 
which is better than what we would get with the standard 
semi-coherent searches [B]. 



V. THE RELATION WITH THE J" STATISTIC 



From (4.151, we see that if we use all SFT pairs avail- 
able, the amplitude sensitivity of the cross-correlation 

— 1/2 

search is proportional to T^^J which is what we would 
get for a fully coherent search. There must thus be a 
close relation between the cross-correlation and the co- 
herent matched filter, and in this section we show that 
this is indeed the case. 

A convenient implementation of the matched filter 
statistic for periodic waves is provided by the so-called 
jF-statistic first defined in for the single interferome- 
ter case, and later generalized to the multi-interferometer 
case in |16j , and a detailed study of the parameter space 
resolution was presented in [33] . Let us start with the 
single interferometer case. 

For defining the J^-statistic, it is convenient to rewrite 
the waveform of (2.2 1. We first separate out the initial 



phase <i>o from the total phase as, 

^(t) = $0 + </?(i) . 



(5.1) 



We decompose the total waveform h{t) in terms of four 
quadratures as. 



1=1 



(5.2) 



where the four amplitudes {^1''} (not to be confused with 
A+ and ^x) are time independent and the {h^} are 



hl{t) : 



a{t) cos Lp{t) 
: a{t) sin (p(t) 



h2{t) 
hi{t) 



b{t) cos (p{t) 
b{t) sin Lp{t) 



(5.3) 



with a(t) and 6(t) defined as in (2.8 1. What this decom- 
position achieves is a separation of the amplitude param- 
eters {/iQ, t, "07 *J'o} from the Doppler parameters. The 
only signal parameters in the quadratures {ft.^} are the 
Doppler parameters while the amplitudes {.4^} depend 
only on the amplitude parameters. 

In order to extract the signal h{t) from the noise, the 
optimal search statistic is the likelihood function A de- 
fined by. 



InA = {x\h) - ^{h\h) 



(5.4) 



df 



(5.5) 



where the inner product (-j-) is defined as: 

The quantity InA is essentially the matched filter and 
is precisely what we should use in order to best detect 
the waveform h(t). An explicit search over the amplitude 
parameters {A^} is avoided by noting that InA depends 
quadratically on the {A'^ } . We can thus analytically find 
the maximum likelihood (ML) estimators {^1^} of the 
amplitudes {^1'^ } by solving the set of four coupled linear 
equations: 

51nA , , 

-OA^^.^^r'^ .= 1,...,4. (5.6) 

The jT-statistic is then defined as the log likelihood ratio 
with the values of the amplitudes {-4^} replaced by their 
ML estimators: 



T := lnA|_4^^_j^ . 
Explicitly, can be written as 

4 B\Fa\^ + A\Fb\^ - CiFaF^ 



(5.7) 



where 

Fa 



FbF*) 



AB-C^ 



Fh^ 



A- 



Tot 



C 



a^{t)dt , B - 
a(t)b{t) dt . 



b^{t) dt 



(5.8) 

(5.9a) 
(5.9b) 
(5.9c) 
(5.9d) 
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We need to write the Fa and Ff, still more explicitly; let 
us start with Fa- We break up the integral for Fa into 
sub-intervals defined by the SFTs, and assume as we have 
been doing all along that a{t) is constant over the SFT 
duration: 

pTi+AT/2 
J Jti-AT/2 

pTi+AT/2 

= Va/ / x{t)e-''f'^*Ut . (5.10) 

J Jti-AT/2 

Writing the phase in a Taylor series around the SFT mid- 
time and keeping the linear terms, we get, 

ifit) = ip{Ti) + i2Tifi{t - Tj) , (5.11) 

which leads to. 



Fq = ^ aie 
I 



l-Ti+AT/2 
Jti-AT/2 

= ^a,e-^(^^)e--/^^i,(//), (5-12) 



and likewise for Fi,. 

Now we are ready to look at T again. From (5.8 1 it is 
clear that is quadratic in the data and from (5.12) it 



is clear that we will end up with an expansion like, 

T = Y,^iJ- (5.13) 



In fact, it turns out that ( |5.13 1 is precisely a linear com- 
bination of the ya defined in (4.1). Explicitly, it follows 
from ( |5.12[ ) that: 



\Fa 



ij 



Similar expressions are obtained for and the cross 

term FaF^ +Fi,F* of the statistic. Combining all of the 
results from above, we see that is a detection statistic 



of the form (4.4) with weights 



uij cx {Abibj + Bajaj - C(a/6j -I- ajbi)) e'^*" . 

(5.15) 

In the case where A B and C ^ A, B, this is seen to 
be proportional to Qjj averaged over cos l and -0 ( |3.17 1 . 
Thus we see that the cross-correlation statistic p is in- 
deed roughly equivalent to the J^-statistic. In principle, 
p using the full signal cross-correlation function Qa from 
( 3. 10 1 , is a function of the Doppler parameters and also of 



{A^,Ax , V'}; this is more like the likelihood-ratio (mod- 
ulo the dependence on the initial phase $0) before max- 
imizing it over the amplitude parameters to obtain the 
J^-statistic. The p calculated with (^Q)cosi,,V' closer to 
the matched filter statistic marginalized over cos t and ijj. 



VI. ESTIMATING THE AMPLITUDE 
PARAMETERS 

Thus far, we have focused on constructing the cross- 
correlation statistic which is optimal for the detecting the 
presence of periodic GWs. Thus, the choice of weights 
given in (4.10) is tailored for measurements of excess 



cross-correlation power, and is not actually an estima- 
tor for the signal amplitude. Estimating the Doppler 
parameters {/o, /i, . . . , n} is easy since we are searching 
over these parameters explicitly. N ote al so that the sig- 
nal cross-correlation function Qa of (3.10) is a function of 
cos L and ip- We could thus, in principle, find the values 
of cos t and V' which maximize p, thus yielding estima- 
tors of these quantities. In practice however, we expect 
it to be more convenient to use a single statistic, such as 
that associated with the averaged Ga given in (3.17), and 



then estimate {A+, A^ , V'} in a follow-up stage.® In this 
section, we show that it is indeed possible to estimate 
{A^, Ax The method presented here is a straight- 
forward generalization of [J? (see also [Sllin]) developed 
for the standard semi-coherent searches. 

The basic idea is to note that the two polarizations 
/i+ and hx appear in the detector with different ampli- 
tude modulations. Therefore, given sufficient measure- 
ments of the it should be possible to extract the 
signal components with different amplitude modulation 
patterns thereby estimating the amplitudes and Ax- 
Let us start by defining the signal cross-correlation func- 
tions Ga ^^^d Ga for two polarizations which are anal- 
ogous to Ga- 



(6.1a) 
(6.1b) 



yi.i - 46 i'i+-t'j+ , 

t'/j — -^e rjx-tjx ■ 

These functions are significant because, just as in ( |4.2[ ), 
they tell us about the mean pa of ya for the two inde- 
pendent polarizations. The contributions of and Ax 
to the mean are respectively: 

(6.2) 



f^l^AlGi and p^^AlG^ 



An estimator of A^ is obtained by minimizing the fol- 
lowing x^-statistic as a function of A"^, 



X 



E 



- AIG + 



The solution of dx^/dA'^ = is easily seen to be, 

-1 

y*aG+ + yc.G+* 



'/3 



E 



(6.3) 



(6.4) 



* Note that the cross-correlations 3^^ are independent of the initial 
phase #0- Thus, it is not possible to estimate $0 if wc restrict 
ourselves to measurements of ya ■ 
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Similarly, the estimator for Ax is, 



X 12 
2 



E 



Since {^^ } depend on the polarization angle -0 through 
the beam pattern functions, both (6.4 1 and ( |6.5[ ) imply a 
search over V'- We expect these estimators to be better 
than the ones used in the standard semi-coherent meth- 
ods simply because it uses a larger number of measure- 
ments including all possible pairs of SFTs. Note that 
these estimators for and are proportional to the 
optimal excess-power statistic p of (4.11), with Qa re- 
placed by ^+ and . 

Finally, while we do not discuss it here, following ^5^, 
this discussion can be generalized to construct a joint 
statistic for A^, and for a general elliptically 
polarized signal. 



VII. PARAMETER SPACE RESOLUTION 

In this section we discuss the parameter space resolu- 
tion required for the cross-correlation statistic p. This af- 
fects the astrophysical significance of the search in terms 
of parameter estimation and also the computational re- 
quirements for carrying out the search. The parameter 
space resolution for a detection statistic p is usually dis- 
cussed in terms of the parameter space metric. This is 
defined as the fractional loss in the signal-to-noise ratio 
when p is calculated at a point in parameter space which 
is slightly different from the point corresponding to the 
actual signal parameters [351 [55] [37] . In our case, we are 
in principle free to consider any subset of all the possi- 
ble SFT pairs in calculating the final detection statistic 
p. However, without some control on which SFT pairs 
are chosen, it seems very hard to get a handle on the 
parameter space metric for the general cross-correlation 
statistic p defined by (4.4). Our suggestion is the follow- 
ing: Choose a time duration T^ax and include only those 
SFT pairs {/, J} for which \Ti-Tj\ < r,„ax. Thus, r,„ax 
can be viewed as the maximum duration over which we 
choose to maintain strict phase coherence. 

If Tinax = ?obs) then we are including all possible pairs, 
and at the other extreme, if Tmax — then we are includ- 
ing only self-correlations and time-coincident correlations 
between different detectors. In the intermediate regime 
the cross-correlation search is closest in spirit to a semi- 
coherent hierarchical scheme which consists of breaking 
up the total data available (say from t = to t = Tobs into 
shorter segments [0, Tmax] , [Tmax, ZTmax] • ■ •• One then 
performs a coherent analysis on each of the segments (us- 
ing, say, the .F-statistic) and combines the results semi- 
coherently [9] [11] [12] . The pair selection criteria would 
lead us to choose all possible SFT pairs within each of the 
segments. Since we have already seen in Sec. [V] that this 
is essentially equivalent to the J^-statistic, the similarities 



between the two schemes is obvious. The two schemes 
are however not exactly identical because this SFT pair 
selection criteria also includes choosing pairs lying in ad- 
jacent data segments (assuming the segments are suffi- 
ciently close to each other). Thus, the cross-correlation 
search with coherence time T^ax will be more sensitive 
than the semi-coherent search with coherent segments of 
duration Tmax but the precise improvement depends on 
the duty cycle of the detectors, i.e. on the gaps between 
the SFTs and the coherent segments. 

With this criteria of choosing pairs, we will see that 
the resolution depends on Tmax the SFT baseline AT. To 
make our results concrete, we will focus on the ground 
based interferometers by taking the frequency range to 
be from 50 Hz to 1000 Hz. Given the similarities with the 
semi-coherent and hierarchical schemes discussed above, 
it is clear that a proper discussion of the metric requires 
a calculation of the parameter space metric for semi- 
coherent searches. This is a combination of the coherent 
metric worked out in detail in [33] [3S], and the semi- 
coherent metric obtained by summing .?^-statistic seg- 
ments. Preliminary calculations have been worked out in 
[S], but a detailed study of its properties is still lacking. 
We will instead resort to order of magnitude estimates 
(which, in spite of their approximate nature, have actu- 
ally turned out to be fairly useful for previous searches; 
see e.g. [T^). 

We can either use the amplitude modulation of the de- 
tection statistic p = Pa by which we mean the varia- 
tion of Pa with a, or we can use the frequency modulation 
reflected in the different frequency bins k and k' used to 
calculate the cross correlation = xljXk'j- Starting 
with the sky-resolution, we identify three factors which 
could be relevant: the detector beam pattern functions, 
the detector-pair baseline Ar/j, and the Doppler infor- 
mation over a duration AT and Tmaxi we discuss all of 
these in turn. The relative importance of these three 
factors depends on the search parameters. 

i. The expectation value of the cross-correlation statis- 
tic varies with the SFT pair index a, and part 
of this variation is due to the geometrical factor 
i/aj -|- 6/6,7 in (3.11 1. Since this variation depends 



on the sky-position, it can in principle be used to 
get sky-position information. The resolution thus 
obtained is roughly given by the angular scales over 
which the beam pattern functions vary. Note that 
this amplitude modulation is due to the rotation 
of Earth around its axis; this is independent of the 
signal frequency and gets mostly averaged out if 
AT is comparable or larger than a day. 

The other reason for the variation of the SNR with a 
is the Afjj term in (3.11). In the case when the 
two SFTs are coincident in time (T/ = Tj), then 
Ar/j is the separation between the two detectors; 
for the LIGO Hanford and Livingston observato- 
ries, this corresponds to a light travel time of about 
10 ms. More generally, the magnitude of Ar/j is the 
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distance between the positions of the two (distinct 
or same) detectors at different times; it could be 
as much as 2AU if Tj — Tj ^ 6 months. On the 
other extreme, it could be zero if we are correlat- 
ing the data with itself (which is what the stan- 
dard semi-coherent methods do); this effect then 
becomes completely irrelevant. If Agw is the wave- 
length of the wave we are trying to detect, the sky- 
resolution associated with Ar/j is inversely propor- 
tional to the frequency: 



A. 



|Ar1 f-\Ar\/c 



(7.1) 



For the Hanford-Livingston pair, this corresponds 
to about 0(60°) at 100 Hz and about 6° at 1000 Hz. 

The third way of getting sky-position information is 
through the Doppler shift. This is only useful if the 
frequency resolution of the SFTs is small enough; 
the maximum Doppler shift is f\v\/c, so for the 
Doppler shift to be important, we must have. 



AT > 



A. 



(7.2) 



The magnitude of Earth's orbital velocity in its or- 



bit is - lO^'^c, so ( 7.2 1 leads to AT > 200 s at 50 Hz 



and AT > 6.67s at 1500 Hz. One relevant baseline 
in this case is the distance traveled by the detector 
in the duration AT. Thus, the sky resolution is 
(see [12] for further details): 



A 



doppler 



gw 



\v\AT 



(7.3) 



For 1800 s SFTs, this corresponds to 6° at 50 Hz 
and 0.2° at 1500 Hz. There is finally the baseline 
corresponding to Tmax, i-e. the distance dmax trav- 
eled by the detector during T^ax- This leads to 



(A0)doppi C 



A„ 



(7.4) 



More generally, the resolution corresponding to 
Tmax (for sufficiently large AT) is precisely the co- 
herent metric calculated in ^33j |38j . 

We see that the first two items above can be viewed as us- 
ing the amplitude modulation information (dependence 
of the SNR on the pair index a) , while the third uses the 
frequency modulation. 

Let us now discuss the resolution in spindown param- 
eters fk- The spindown term in A<i>/j appears in the 
combination fk{T^^^ - T^+^). Thus, it is clear that for 
T/ 7^ Tj this leads to a spindown resolution of. 



1 



maxj,4|T|-+i-T)+i|} 



(7.5) 



Thus, if we were to consider all possible pairs from a 
given set of SFTs, and if we define the reference time to 



be in the mid-point of the observation duration, then we 



-(fc+i) 

bs 



would have 5fk oc T^ 

We can also consider the frequency resolution {6f)sft 
(AT)^^ of the SFTs themselves. The corresponding res- 
olution in fk is defined by the smallest change in fk re- 
quired to change the frequency by a (<5/)sft over the full 
observation time Tobs- This leads to 5fk = {^f)sh/T^bg 
for fc= 1,2.... 

Let us conclude this section by giving a short numerical 
example for the case when we correlate data from a pair 
of spatially separated detectors at the same times. We 
consider frequencies of 100 Hz and 1000 Hz, and two sky 
positions: one at the celestial equator and one at 45° de- 
grees above it. In each case we consider sources with the 
optimal orientation t = 0, without any spindown parame- 
ters, and with -0 = 0. The total observation time is taken 
to be Tobs = 1 yr and the SFT basehne is AT = 30 min. 
We assume that the two data streams are coming from 
the LIGO Livingston and Hanford interferometers. For 
performing the cross-correlations, we use. 



Q{t;f,n) = X{t;n){g{n)) cos L,-4,, 



(7.6) 



where, X(t; n) is a proportionality constant. We con- 
sider essentially identical time segments - same barycen- 
tric time - in the two detectors. In a year's worth of 
observation time there are little over 17,000 such time 
segments, each of 30 minutes duration. Thus the time- 
segment indices /, J each, sequentially run over the full 
observation time. The relevant quantities Q, Q and A in 
(7.6) which carry the same indices also do the same over 



the observation time - thus we may think of each of them 
as functions oi t - the segment time-stamp; thus / or J 
is replaced by t. 

For the signal only case, the cross-correlation can be 
written explicitly as: 

B{n,n') ^ K{n) I (^(t; n))cost,^/i(i)(i; «')^(2)(i; "') 

Jo 

(7.7) 

where the subscripts in and h(2) refer to the two 
distinct detectors we are considering. We have chosen. 



A-i(n) = 



1 

AT 



(7.8) 



We choose the proportionality constant A such that it 
is inversely proportional to square of the average total 
power accessible to the network for a particular direction 
of the sky in the interval AT of the SFTs. Thus we have. 



\-\t-n) = AT{g{t-n)) 



2 

cos L^tp' 



(7.9) 



This is in the spirit of the normalization scheme adopted 
in [T5]. Figure [l] shows B{n,n') evaluated numerically 
for point sources at different positions. We note that 
the maximum value of B is 5. This is the result of the 
average value of Q we have chosen in defining the filter 
function together with the fact that we have chosen op- 
timally oriented sources for the numerical computation. 
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The sky-resolution is characterized quantitatively by the 
FWHM (full width at half maximum) of the PSF. From 
the figures it turns out to be ~ 8° for /q = 1000 Hz and 
~ 80° for /o — 100 Hz. We observe that the agreement 
between the order of magnitude estimates obtained ear- 
lier and the actual values computed from the figure is 
satisfactory. 



VIII. DISCUSSION 

We summarize the main results of this paper. We have 
generalized the cross-correlation statistic, traditionally 
used for the stochastic gravitational wave background 
searches, to periodic gravitational waves. The features of 
periodic waves, not present in the stochastic background 
signals, are non-stationarity and long-term coherence. 
The non-stationarity may need to be taken into account 
depending on the frequency resolution, and the long-term 
coherence implies that we can in principle cross-correlate 
data segments from arbitrary times and arbitrary detec- 
tors. This makes the method very flexible, and these are 
some of the possibilities: 

i. We can, if we wish, correlate all possible short data 

segments. If this is done, then we showed that 
the resulting detection statistic is very close to the 
.F-statistic corresponding to a full matched filter 
statistic. This is ideally the most sensitive method, 
but it's computational cost becomes prohibitive for 
wide parameter space searches. 

ii. At the other extreme, we can choose to correlate only 

data segments taken from distinct detectors at the 
same (or very close) times. This is the closest 
in spirit to the standard directed stochastic back- 
ground searches using aperture synthesis. In this 
mode of operation, the search is not computation- 
ally intensive, and is very robust against signal un- 
certainties. However, this also implies poor resolu- 
tion in parameter space, and thus more expensive 
follow-ups to verify possible detections and to esti- 
mate the signal parameters. 

iii. From the perspective of this paper, the standard 

semi-coherent searches such as PowerFlux, Stack- 
Slide and Hough all correspond to the special case 
in which we consider only self-correlations. The 
procedure of considering weighted sums of the 
cross-correlation power is closest to the PowerFlux 
method. In fact, many of the lessons learnt in the 
PowerFlux searches should be applicable here with 
suitable modifications. For example, the estima- 
tion of the signal amplitudes developed originally 
for PowerFlux carries over rather straightforwardly. 

iv. In intermediate regimes when we correlate data seg- 

ments separated by a maximum coherence time 
Tmax < Tobs, the cross-correlation search is similar 



to a hierarchical search in which we combine seg- 
ments of demodulated data. Though, as discussed 
in Sec. |VII[ there are differences between the two 
with the cross-correlation search being somewhat 
more sensitive. 



Conceptually, this method thus provides a unified frame- 
work for all the known periodic wave searches, and this 
might be useful in various calculations and applications. 
Each of the above modes of operation correspond to tun- 
ing the maximum coherence time all the way from small 
values to the total observation time. The precise value 
chosen for a specific appHcation depends on the trade-offs 
between computational cost, sensitivity, and robustness 
against signal uncertainties. The additional parameter 
which figures importantly in this trade-off is the length 
AT of the short data segments. 

There are a number of open issues for future work. 
An important question is to get a detailed understand- 
ing of the trade-offs mentioned above for various types 
of searches including all sky searches for isolated GW 
pulsars, signals from known binary systems or from in- 
teresting areas such as the galactic center etc. This will 
help us better decide how to best use our computational 
resources and to maximize our chances of making a de- 
tection. Another important issue, which feeds into this 
optimization problem, is to study the general parameter 
space metric. To date we only have a proper understand- 
ing of the coherent metric, i.e. case (i) above. For the 
other cases, we have estimates of the parameter space 
resolution and which are often sufficient for many appli- 
cations, but a full understanding is still lacking. In ad- 
dition, it would be interesting to compare the estimation 
of the amplitude parameters Ax} (and ^) obtained 
from (6.4 1 and (6.5) with the maximum likelihood esti- 
mators obtained from the JF-statistic calculation. In the 



limit when we consider all possible correlations, we would 
expect the two estimates to be very close to each other. 
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(a) Source ata = , 6 = ; f =1000 Hz 



(b) Source at a = 0°, 6 = 45°: fg=iooo Hz 





(c) Source ata = , 6 = ; t = 100 Hz 



(d) Source at a = , 6 = 45 ; f =100 Hz 





FIG. 1: The point spread functions (PSFs) for sources with frequencies 100 Hz [(c) and (d)] and 1000 Hz [(a) and (b)]. The 
source is taken at the celestial equator [(a) and (c)] and 45° above the celestial equator [(b) and (d)]. In all the cases, the 
source orientation is taken to be optimal, i.e., (, = '0 = 0. 



APPENDIX A: INCLUDING 
SELF-CORRELATIONS AND 0{hl) 
CORRECTIONS 

In this section we relax the two assumptions of only 
looking at for I ^ J and h <^ n. We allow self corre- 
lations (which, by themselves, are used in the standard 
semi-coherent searches), and we keep terms oiO{h^) but 
still neglect C(/iq) terms. 

Let us again start from the general statistic p defined 
in (4.4 1, and let us calculate its mean and standard devi- 
ation with the two assumptions relaxed. In general, we 
have = y,ji so that yjj is real and so is the corre- 
sponding weight uji; yjj is in fact just the power in a 
single SFT bin. We will denote yjj simply by 3^/ and uji 
by uj. 

The mean is easy to calculate: 



1 



fj-ij 



2AT 



SiSij + hlGu 



(Al) 



Thus, the mean is non-zero in the absence of a signal only 
for the self-correlation terms. In general, p will contain 
self-correlations, and also correlations of distinct pairs. 
However, we want to be completely general and we do 
not assume that it contains all the possible pairs. This 
is then the expression for the mean: 



{p) := /i = 



AT ^ 



uiSi + hlY^iuJ^+ulGl) . (A2) 



It is to be understood that the first sum in this equation 
only contains the self-correlations and the second sum 
contains all the SFT pairs we have chosen to include, 
including the self-correlations. 



The variance calculation is somewhat more involved. 
Before looking at the variance of p itself, let us look at 
{yijyKL)- Note that for the pure noise terms: 



(A3) 



Here, we use the notation that indices within parentheses 
are symmetrized over: ^(/j) = {Xjj + Xjj)/2. This also 
covers the I = J = K = L case, so there is no need to 
consider that separately. 

Consider now the signal. In general, the terms in 
{yijyKh) with odd powers of h will vanish because the 
noise is assumed to have zero mean. Thus, schematically, 
we will have 



{yuynL) = A + Bhi + chi 



(A4) 



Let us ignore the /iq terms and focus only on the second 
order terms. The reader can convince herself that we 
only need to keep the following terms in yijyKL'- 

(A5) 



Putting together ( A3 1 and ( A5 1 , we end up with 



"o 
AT 



GiijSDKSi''^ +Sj(^jgL^KSi'A . (A6) 



We are now ready to look at the variance of p. Let us 
define pa — Uq3^q + u* so that p = J2a Pa- Then, we 
have 



Var (p) 



|Var(p„) 



E 



Gov {pc^^pp) . (A7) 
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Let us start with the variances 

Var (pij) = (pIj) - 



(A8) 



For I ^ J, pij — 0{hl) so that p.jj can be ignored. 
Thus, in this case we get: 



Var(p/j) = 2|u/j|^ 



oil) c,(J) 
On On 

4Ar2 



For the I = J case, we can no longer ignore the /io, term. 
Keeping terms up to O^Hq) we end up with 



^1 



Var {pj 




(AlO) 

Turning now to the covariances, first note that if 
I,J,K,L are all distinct, then up to ©(/iq) terms, 
Cov {pu, pkl) — 0; thus we need at least one pa ir of 
matching indices to get a non-zero result. Using (A6) 



the expressions for all the non-zero cases are the follow- 
ing (/ 7^ J) ignoring, as always, the 0{hQ) terms: 



{yiiyij) = ^{gij + gji)s^j^ , 



(Alia) 
(AUb) 



(yiiyjj) 



oil) qiJ) 

On On 

4AT2 



2AT 



(Allc) 



It turns out that the only non-zero covariance is 

Gov {pi,pij) = -^uiSi'^ (uuGh + uljGij) . (A12) 

We are almost done now. Substituting the results of 
(|A9l), (|ATo|, and (|A12|) in (|A7| we get 



E \nij\'{gjsi;'^ + gjsi^^ 

*ijgij)si' 



(A13) 



Here we have defined the variances in the absence of a 
signal: 



'■(o)J 



f(0),/J 



2Ar2 

On On 

4AT2 



(A14a) 



(AUb) 



It is convenient to write (A13l in the abbreviated form 

(A15) 



2 2 I 1,2 2 



'(0) 



'0"(1) 



where the definitions of CTq ^J^d cr^i) are obvious from 
dATl. 



We are finally ready to derive the equation for the sen- 



sitivity, i.e. the analogs of (4.9 1 and (4.151. (4.7 1 for 



the threshold is unchanged as long as we use (T(o) instead 
of (7 in that equation^. (4.8 1 for the false dismissal rate 
becomes: 



Keeping terms linear in /iq, we get 



erfc-i(27) = ^^^-''°"'^^ 



(A16) 



/2cr 



(0) 



= erfc-i(2a) - Y.^uJ^ + u*Jl) 

V2cr(o) „ 



erfc-\2a) . (A17) 



Solving for ho leads to the generalization of (4.9 1 



Ea(""^" + "a^a) + ^fi) Crfc" ^ (2a) /%/2cr(o) 

(A18) 

Finding the optimal weights is now not as straightforward 
as before. However, we note that when cr(]^-) is ignored, 
then the optimal weights are again given by (4.10) ex- 



cept that now it holds also for the self-correlations. In 
the general case when we do not ignore <J(i), it is simpler 
to continue using the optimal weights derived earlier in 
(4.101, and to substitute it in (AI8I to derive the corre- 



sponding sensitivity. 
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